Hydrologic Design 4501

Chapter 5: Infiltration

Ardeshir Ebtehaj
University of Minnesota
Table of Contents
Chapter 5: Infiltration 1- Soil Properties 2- Vertical Distribution of Subsurface Water 3- Saturated Flow in Porous Media 3-1 Hydraulic Head 3-2 Darcy’s Law for Saturated Flow 4- Unsatuarated Flow in Porous Media 4-1 Darcy’s Law for Unsaturated Flow 4-2 Soil Water Characteristic Curve 4-3 Richards' Equation 5- Infiltration Models 5-1 Conceptual Understanding of Infiltration 5-2 Philip's Model 5-3 Green-Ampt Model

1- Soil Properties

Soil is a porous medium consisting of solid grains and interconnected voids that can be filled with air or water.
The ability of soil to retain and transport water depends on the distribution of its particle size, pore openings and the chemical composition of the soil particles. Particle Size Distribution (PSD) is a common method to classify soils based on some standard criteria.
For large grains (D > 0.05 mm), a series of successively smaller screens is used to sieve out the mass percentage above each screen size. For small or fine grains (D ≤ 0.05 mm), the settlement rate in a liquid is used to determine the PSD.
Figure 1: Typical sieves used in creating a PSD (left). A typical hydrometer test for fine grain PSD analysis (middle). A sample PSD for three different soils (right).
The primary particle size classification used by hydrologists is from the United States Department of Agriculture (USDA) (shown below). The range of grain sizes in a soil determines the soil texture, which is also defined by the USDA soil texture triangle, as shown below.
For example, a soil with loam texture is defined as 40% sand, 40% silt, and 20% clay. It is important to note that, fine grain soils (silt and clay) typically have more total pore spaces due to their larger surface area and mineral composition. However, as we will find out, this does not necessarily result in a larger water mass transport.
Figure 2: USDA soil particle (left) and soil texture (right) classification.
From a hydrologic point of view, we are interested in the water holding capacity and water transportation in soil. Therefore, we need to define some basic mass and volume relationships to specify the amount of air, water and soil particles in a control volume.
Bulk Volume (V): The total control volume of the soil matrix is , where is the volume of voids and is the volume of solids.
Volume of Voids (): The total pore space of the soil matrix, , where is the volume of air and is the volume of water.
Figure 3: Classifying soil as a three phase system (air, water and solids) to create mass and volume relationships (Credit: the COMET Program and Hillel 1998).
Porosity (η): Dimensionless measure of void space in soil matrix:
,
where is the mass of the solids, is the bulk density, which typically ranges from 1.0 to 1.6 [], and is the particle density, which typically ranges from 2.6 to 2.75 [].
Note: There is a parameter called effective porosity (), which excludes isolated pores and pore volume occupied by capillary water adsorbed on clay minerals or other grains. The effective porosity is typically less than total porosity and better explains the soil water that contributes in water transport in the soil matrix. We should emphasize that the total porosity is the total void space in the soil whether or not it contributes to fluid flow.
Soil Volumetric Water Content (θ): Dimensionless measure of the soil moisture content,
where .
Relative Saturation (S): Volumetric water content normalized by porosity,
where .
Figure 4: Mean (standard deviation) of porosity as a function of soil type (left, Clapp & Hornberger 1978). Mean (mean±one standard deviation) of porosity and effective porosity as a function of soil type (right, Rawls et al. 1983).

2- Vertical Distribution of Subsurface Water

In hydrology, it is useful to conceptualize different zones in the subsurface soil, based on typical soil moisture profile as shown in Figure 5. The saturated zone refers to sections where pores are completely filled with water and the pressure is above atmospheric pressure. This zone is typically what we think of as groundwater and follows simpler saturated flow equations. The groundwater table is where pressure in soil is equal to the atmospheric pressure.
The vadose or unsaturated zone is the partially saturated soil directly above the groundwater table. The pore pressure in this region is less than atmospheric pressure. This zone is governed by more complex flow equations than the saturated zones and can be split into three additional sub-zones.
Figure 5: Classification of vertical soil structure based on the moisture content
Figure 6: Capillary rise in soils with different textures after 72 days (Lohman, 1972)

3- Saturated Flow in Porous Media

3-1 Hydraulic Head

Water flows down the prevailing energy gradient, which we can be defined at a given point in a fluid flow field in terms of:
Potential Energy: Kinetic Energy: Pressure Energy: .
Thus, the total energy of the fluid flow is as follows:
[Joules].
However, it is often convenient to express the above components in terms of hydraulic head, which is total energy per unit weight of a water parcel as follows:
[].
Here, the pressure head (P) is the gauge pressure (e.g., ), which means that we should not account for atmospheric pressure.
In subsurface flow, the velocity is typically very slow and we can assume , therefore we have
[].
Therefore, for saturated flows, the hydraulic head driving the flow is dependent on gravitational component (z) and hydrostatic pressure forces ().

3-2 Darcy’s Law for Saturated Flow

In 1856, Henry Darcy ran experiments similar to Figure 7 to test sand filters for purifying the drinking water of Dijon in France.
Figure 7: Left: Darcy's experimental set up ()
Darcy found that the steady state volumetirc flux of water discharge per unit area through a saturated porous medium is proportional to the hydraulic gradient across the filter length as follows:
[],
where [] or is a constant of proportionality known as the saturated hydraulic conductivity.
As an example, in Figure 7, we can expand the Darcy's law as follows:
.
However, and thus,
.
As shown in Figure 7, the soil hydraulic conductivity linearly transforms the hydraulic gradient to a flux depending on the soil texture. Some typical values are shown in the following table, which vary by orders of magnitude. As a general rule, soils with higher clay content have much lower values, because of a greater piezometric head loss due to more tortuous flow path through the smaller pores.
Table 1: The hydraulic conductivity for different USDA soil textures (Clapp & Hornberger, 1978)

4- Unsatuarated Flow in Porous Media

4-1 Darcy’s Law for Unsaturated Flow

Darcy’s Law was originally created for saturated flow below the water table, however, many important hydrologic processes occur in unsaturated vadose zone.
In this unsaturated zone, we must deal with pressure that is often below atmospheric pressure due to the nagative suction or capillary pressure .
Recall from fluid mechanics, the attraction of water to a surface (adhesion) and to itself (cohesion) causes water to be pulled up a narrow capillary tube until that capillary force is balanced by the downward gravitational force on the water column.
The height of capillary rise is represented by:
where σ is the water surface tension .
Figure 8: Conceptual (left) and actual capillary tubes (right).
We can think of the tortuous pore spaces in soils as a series of twisted capillary tubes with negative pressure. The smaller the pores, the higher the capillary rise. Therefore clayey soils will have a larger capillary fringe height than sandy soils (Table 2).
Table 2: Capillary rise in soils with similar porosity after 72 days (from Lohman, 1972)
Let us represent the capillary head as a negative pressure head:
[]
Thus, we can modify what we consider as the hydraulic or piezometric head as follows:
.
Now let us define the Darcy’s Law for unsaturated flow with this new definition of piezometric head. For simplicity, we are going to define it only in z-direction that is our primary focus for studying infiltration.
We know from experience that the unsaturated hydraulic conductivity and capillary head are functions of soil moisture content θ.
Thus, we can use the chain rule to rewrite the Darcy’s Law as follows:
Now for simplicity, we can define the first term on the left hand side as a function of θ called the the soil water diffusivity with unit to get our final form of Darcy’s Law for unsaturated flow:
[]
It is important to note that in the above representation, the Darcy’s law is expressed as a function of soil moisture gradient, which is easier to measure than the soil capillary head.
The moisture flux in the unsaturated zone is driven by two factors:
  1. Metric forces from wetter to drier soils.
  2. Gravity, whieh is always downward.
To get a better understanding of the unsaturated flow, the figure below shows the signs of these terms in the above equation for two ideal scenarios. Note that, here we define z as positive upward.
Figure 9: Conceptual sketch of dominant forces in unsaturated moisture flow in a vertical soil column.

4-2 Soil Water Characteristic Curve

In order to use the unsaturated Darcy Law, we need to quantify functional dependence of both and to soil moisture often from lab tests. Note that here after we drop the subscript z from for notational convenience.
The soil water characteristic curve (SWC) or retention curve is the relationship between and θ.
Figure 10 shows that the capillary head increase as the soil moisture decreases. The unsaturated hydraulic conductivity as a function of soil moisture is determined similar to by measuring outflow for as set of capillary head differences and unsaturated soil moisture content. Note that the capilalry head approaches to zero as the soil becomes saturated and tends to .
Figure 10: Schematic representation of changes in capillary pressure (absolute values) and hydraulic conductivity as a function of relative saturation . When the soil becomes saturated, the hydraulic conductivity approaches to the saturated hydraulic conductivity .
There have been multiple attempts to parameterize the SWC and a function of soil texture and moisture. One of the widely used methods in hydrology is the Brooks-Corey model as follows:
Table 3: Saturated matric head and Brooks-Corey parameter (b) (from Clapp & Hornberger, 1978). The values in parentheses represent the uncertainty in terms of one standard deviation.
Now based on the provided discussion, we can have a more precise explanaiton of soil moisture field capacity and wilting points.
Field Capacity : This is the maximum water content that can be retained by a soil against gravitational forces. A practical explanation of this quantity is the moisture content of a soil 2 to 3 days after rainfall or irrigation. The field capacity is formally defined as the soil water content corresponding to the [] capillary pressure, which can be obtained using the Brooks-Corey model.
Wilting Point : This is the defined as the water content at which most plants will begin to wilt since they cannot access the water they need from the soil. It is formally defined as the soil water content at [].
Available water content : The difference between field capacity and the wilting point, which defines how much water plants have access to uptake from the soil:
.
Figure 11: Left: Illustration of various wetness conditions. Right: Graphical representation of wetness conditions for different soil textures (Credit: The COMET Program)

4-3 Richards' Equation

Thus far, we have been discussing steady state flow equations, however, in order to define the governing equation of unsaturated flow, we need to incorporate the time through the conservation of mass. We can derive this by looking at a control volume and assume that the only flux is in the z-direction:
Figure 12: Control volume element used to derive continuity equation for unsaturated flow.
Accumulation of (mass/time) =
Accumulation of (mass/time) =
Input (mass/time) =
Output (mass/time) =
Source/Sink =
Note that all the terms have units of mass per time such as []. [] is a term that can be used to represent plant water uptake as a sink .
Combining these equations leads to the following continuity equation:
.
Note: This conservation equation has unit of [] since it is .
Finally, we are able to substitute the unsaturated Darcy’s Law (momentum equation) for the in the above continuity equation to obtain the one-dimensional Richard’s Equation that explaines the chnages of moisture dynamics throughout the soil column,
[].
Note: , , and are functions of z and thus we cannot take out of the derivative operator. This makes the Richards’ equation a non-linear equation. Since is extremely nonlinear with respect to θ. The solution of the Richards' equation through numerical methods is often challenging and subject to numerical instabilities.

5- Infiltration Models

We have now derived a detailed equation for unsteady, unsaturated flow in porous medium. However, in hydrology, we are typically interested to estimate how much precipitation water infiltrates at the surface into the vadose zone. Given the difficulties in numerical solution of the Richards’ equation over a large watershed, it is more practical to make simplified assumptions for calculating the infiltration fluxes.

5-1 Conceptual Understanding of Infiltration

From earlier discussion, the volumetric infiltration flux at the earth’s surface can be written as follows:
[ or ]
where is shorthand for at and time t. To reiterate, has a negative sign when z is positive upward and has a positive sign when z is positive downward.
We can then simply define the cumulative infiltration by integrating the infiltration rate over time as follows:
Based on this understanding, let’s define three infiltration scenarios that may occur during a precipitation event of rainfall rate p [].
For brevity, we will omit the subscripts nowing that we discussing the downward water flux at and .
1) Ponding after a certain time ():
where is the time to ponding after which water accumulates on the ground and the soil at is assumed to be saturated. Note that after the ponding time, we assume the values of and correspond to the saturated values. As the time passes, the diffusive element, diminishes leaving only the gravitational component . In other words, at the surface as .
2) Immediate Ponding ():
3) No ponding ():
Figure 13: Typical water content profile (left) and evolution of infiltration rate (right)

5-2 Philip's Model

To approximate the infiltration rate, Philip (1960) provided an approximate solution to Richard’s equation as follows:
,
Where is called the sorptivity coefficient:
[]
where is the capillary head at the edge of the wetting front defined by:
and is the saturated capillary head, b denotes the Brooks-Corey’s constant.
Integrating the Philip’s equation yields the cumulative infiltration function as follows:
.
Figure 14: Phillip’s model assumes that infiltration rate decays as square root of t. we need to note that the Philip’s model is undefined at as .
----------------------------------------------------
Example Probelm 5.1: A small horizontal tube with a cross-sectional area of 40 is filled with soil. The front end of the tube is saturated, and after 15 minutes 100 of water have infiltrated into the cube. If the saturated hydraulic conductivity of soil is 0.4 , determine how much infiltration would have taken place in 30 minutes if the soil column had initially been placed upright with its upper surface saturated.
Solution: The cumulative infiltration depth in the soil column is . In a horizontal infiltration, cumulative infiltration is only a function of soil suction head. Therefore, after , we have
For infiltration down a vertical column we have
----------------------------------------------------

5-3 Green-Ampt Model

Figure 15: Conceptual diagram of the wetting front and variable used in the Green-Ampt model.
For the Green-Ampt model, we assume a sharp wetting front (Figure 15) and we begin with Darcy’s Law for saturated flow:
where is the length of the wetting front and H is the ponding depth at the surface and is the infiltration rate.
Now using this equation, we can make a few simplifying assumptions in order to calculate the time to ponding if the rainfall rate p is less than infiltration capacity at (no ponding).
Let’s substitute p for and assume that H is negligible.
(1)
At the same time, we can compute that the cumulative depth of infiltrated water at time t, defined with a sharp wetting front, as follows:
(2)
At the same time, the cumulative infiltration at the ponding time can be represented as .
Therefore, at the ponding time, from equations (1) and (2), one can have
,
and the ponding time can be computed as follows:
Now, let’s develop an equation for cumulative infiltration when . There are two separate equations that can be set equal if we assume H is negligible:
.
Now, we also know that
.
Then, upon integration, we can obtain for :
.
Note that the above equation is implicit and its solution requires an interative rrot find approach.
The cumulative infiltration as a function of time can be used to calculate runoff as follows:
A Summary of the Green-Ampt Model:
1- Ponding time:
[], where and p is the rainfall rate.
2- Cumulative infiltration before ponding time:
.
3- Cumulative infiltration at ponding time:
4- Cumulative infiltration at :
:
5- The infiltration rate at time t can be computed as follows:
----------------------------------------------------
Example Probelm 5.2: Compute the infiltration rate and cumulative infiltration after one hour of infiltration into a sandy loom soil that initially had an effective saturation of 30%. Assumed water is pounded with a small depth on the surface, saturated soil moisture [-], , .
Solution:
and
The cumulative infiltration at is calculated by finding the root of the following equation
we can generally begin with and iterate, which results in . The infiltration rate is
.
clear
% Inputs
theta_s = 0.486; % saturated soil moisture [-]
s = 0.3; % effective initial saturation ratio [-]
psi_f = 16.7; % wetting suction head [cm]
K_sat = 0.65; % Saturated hydraulic conductivity [cm/hr]
t = 1; % time [hr]
% Outputs
D_theta = (1-s)*theta_s;
inf_func = @(F) F - K_sat*t-psi_f*D_theta*log(1+F/(psi_f*D_theta));
% fzero(fun(x),x0) finds the root of fun(x) starting the search at x0
F_est = fzero(inf_func,2*K_sat*t);
disp(['Total infiltration F = ',num2str(F_est),' [cm]'])
Total infiltration F = 3.1672 [cm]
f = K_sat*(psi_f*D_theta/F_est+1);
disp(['The infiltration rate after 1 hour f = ',num2str(f),' [cm/hr]'])
The infiltration rate after 1 hour f = 1.816 [cm/hr]
----------------------------------------------------
Example Probelm 5.3: Compute the ponding time and the depth of water infiltrated at ponding for silt loam soil with 30% initial saturation subject to rainfall intensity of (a) 1 and (b) 5 and then for the latter rainfall rate calculate the cumulative infiltration and its rate after 1 . Assumed water is ponding with a small depth on the surface, saturated soil moisture [-], , .
Solution:
(a) For :
(b) For :
clear
% Inputs
theta_s = 0.486; % saturated soil moisture [-]
s = 0.3; % effective initial saturation ratio [-]
psi_f = 16.7; % saturated suction head [cm]
K_sat = 0.65; % Saturated hydraulic conductivity [cm/hr]
D_theta = (1-s)*theta_s;
% Solution
% (a)
i_r = 1; % rainfall intensity [cm/hr]
t_p = K_sat*psi_f*D_theta/(i_r*(i_r-K_sat)); % [hr]
F_tp = i_r*t_p; % [cm]
% (b)
i_r = 5; % rainfall intensity [cm/hr]
t_p = K_sat*psi_f*D_theta/(i_r*(i_r-K_sat)); % [hr]
F_tp = i_r*t_p; % [cm]
% cumulative infiltration and infiltration rate after t = 1 [hr]
t = 1; % time [hr]
inf_func = @(F) F - F_tp - psi_f*D_theta*log((psi_f*D_theta+F)/(psi_f*D_theta+F_tp))-K_sat*(t-t_p);
F_est = fzero(inf_func,2*K_sat*t);
disp(['Total infiltration F = ',num2str(F_est),' [cm]'])
Total infiltration F = 3.0179 [cm]
f = K_sat*(psi_f*D_theta/F_est+1);
disp(['The infiltration rate after 1 hour f = ',num2str(f),' [cm/hr]'])
The infiltration rate after 1 hour f = 1.8736 [cm/hr]
----------------------------------------------------